%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% COMPUTE FARMERS' LOST SURPLUS FROM ALTERNATIVE POLICIES

% "DEFORESTATION IN THE AMAZON:
% A UNIFIED FRAMEWORK FOR ESTIMATION AND POLICY ANALYSIS"

% by Eduardo Souza-Rodrigues

% This version: November 2018

% OBSERVATIONS: 
% Run this program after running "dem_def_reg.m" & "dem_def_compute_demand.m"

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% (A) PREPARE DATA

% I. Load data

% Load Data Farm (calculated in "dem_def_reg.m")
% Data_farm = [Code Y Area Qa1 Qa2 Uivqr beta_ivqr(1,Uivqr)'];
load Data_farm_1;
load Data_farm_2;
load Data_farm_3;
load Data_farm_4;

if penalty == 0      % Actual Monitoring (Main Specification)       
   
    % Load Total Demand per farm size (calculated in "dem_def_compute_demand.m")
    % Demand_farm = [transfer' tot_demand_iv' tot_demand_ivqr']
    % Total Demand per farm size - Crops & Pasture
    load Demand_farm_1cp;
    load Demand_farm_2cp;
    load Demand_farm_3cp;
    load Demand_farm_4cp;
    
    % Load Demand per farm size & per municipality (calculated in "dem_def_compute_demand.m")
    % deforest_farm_j = [Code Area Demand] 
    % OBS: Demand is a matrix whose lines correspond to municipalities and columns correspond to different taxes
    % Demand per farm size & per municipality - Crops & Pasture - IVQR
    load deforest_farm_1cp;
    load deforest_farm_2cp;
    load deforest_farm_3cp;
    load deforest_farm_4cp;
    
    % Demand per farm size & per municipality - Crops & Pasture - 2SLS
    load deforest_farm_1cp_iv;
    load deforest_farm_2cp_iv;
    load deforest_farm_3cp_iv;
    load deforest_farm_4cp_iv;
    
elseif penalty == 1     % Fines post-2004 at pre-2004 level         
    
    % Restrict Calculations to Models that Include Fines Among Covariates 
    % I.e., exclude model specifications: 
    % model = 3 (no distance to ibama, nor fines), and
    % model = 5 (no fines)
    if model ~= 3 && model ~= 5
        
        % Load Total Demand per farm size (calculated in "dem_def_compute_demand.m")
        % Demand_farm = [transfer' tot_demand_iv' tot_demand_ivqr']
        % Total Demand per farm size - Crops & Pasture
        load Demand_farm_1cp_fines;
        load Demand_farm_2cp_fines;
        load Demand_farm_3cp_fines;
        load Demand_farm_4cp_fines;

        % Load Demand per farm size & per municipality (calculated in "dem_def_compute_demand.m")
        % deforest_farm_j = [Code Area Demand] 
        % OBS: Demand is a matrix whose lines correspond to municipalities and columns correspond to different taxes
        % Demand per farm size & per municipality - Crops & Pasture - IVQR
        load deforest_farm_1cp_fines;
        load deforest_farm_2cp_fines;
        load deforest_farm_3cp_fines;
        load deforest_farm_4cp_fines;

        % Demand per farm size & per municipality - Crops & Pasture - 2SLS
        load deforest_farm_1cp_iv_fines;
        load deforest_farm_2cp_iv_fines;
        load deforest_farm_3cp_iv_fines;
        load deforest_farm_4cp_iv_fines;
        
    end
    
end


% II. Define Variables

% Total Private Area
private_area  = sum(Data_farm_1(:,3)) + sum(Data_farm_2(:,3)) + sum(Data_farm_3(:,3)) + sum(Data_farm_4(:,3));

if penalty == 0      % Actual Monitoring (Main Specification)        
    
    % Transfers 
    transfer_cp  = Demand_farm_1cp(:,1);
    
    % Main Demand for Deforestation - IVQR
    demand_def_cp  = Demand_farm_1cp(:,3) + Demand_farm_2cp(:,3) + Demand_farm_3cp(:,3) + Demand_farm_4cp(:,3);

    % Demand for Deforestation - 2SLS
    demand_def_cp_iv  = Demand_farm_1cp(:,2) + Demand_farm_2cp(:,2) + Demand_farm_3cp(:,2) + Demand_farm_4cp(:,2);

elseif penalty == 1     % Fines post-2004 at pre-2004 level         
    
    % Restrict Calculations to Models that Include Fines Among Covariates 
    if model ~= 3 && model ~= 5

        % Transfers 
        transfer_cp  = Demand_farm_1cp_fines(:,1);

        % Main Demand for Deforestation - IVQR
        demand_def_cp  = Demand_farm_1cp_fines(:,3) + Demand_farm_2cp_fines(:,3) + Demand_farm_3cp_fines(:,3) + Demand_farm_4cp_fines(:,3);

        % Demand for Deforestation - 2SLS
        demand_def_cp_iv  = Demand_farm_1cp_fines(:,2) + Demand_farm_2cp_fines(:,2) + Demand_farm_3cp_fines(:,2) + Demand_farm_4cp_fines(:,2);
        
    end

end    


%% (B) LOST SURPLUS FROM UNIFORM TAX

% I. IVQR

% Close to prop
close = 0.015;

% Find Point at the Demand Curve in which the share of deforest is close to "prop"
index1      = find((demand_def_cp./private_area >= prop - close) & (demand_def_cp./private_area <= prop + close));
[~, index2] = min(abs((demand_def_cp(index1)./private_area) - prop));

% Value of the uniform tax
tax_ivqr = transfer_cp(index1(index2));

% Trapezoid Area Below Demand
area_below_demand_cp = trapz(transfer_cp(1:index1(index2)),demand_def_cp(1:index1(index2)));

% Tax Revenue from Uniform Tax
tax_revenue_cp = transfer_cp(index1(index2))*demand_def_cp(index1(index2));

% Lost Surplus from tax
LS_uniform_tax_cp = area_below_demand_cp - tax_revenue_cp;


% II. 2SLS

% Find Point at the Demand Curve such that share of deforest is close to "prop"
index1      = find((demand_def_cp_iv./private_area >= prop - close) & (demand_def_cp_iv./private_area <= prop + close));
[~, index2] = min(abs((demand_def_cp_iv(index1)./private_area) - prop));

% Value of the uniform tax
tax_iv = transfer_cp(index1(index2));

% Trapezoid Area Below Demand
area_below_demand_cp_iv = trapz(transfer_cp(1:index1(index2)),demand_def_cp_iv(1:index1(index2)));

% Tax Revenue from Uniform Tax
tax_revenue_cp_iv = transfer_cp(index1(index2))*demand_def_cp_iv(index1(index2));

% Lost Surplus from tax
LS_uniform_tax_cp_iv = area_below_demand_cp_iv - tax_revenue_cp_iv;


%% (C) LOST SURPLUS FROM 80% RULE


% I. Define Variables

% Tax/Transfer
tax_cp = transfer_cp';

if penalty == 0
    
    % Sample size for each farm size
    T(1) = size(deforest_farm_1cp,1);
    T(2) = size(deforest_farm_2cp,1);
    T(3) = size(deforest_farm_3cp,1);
    T(4) = size(deforest_farm_4cp,1);
    TT   = max(T);

    % Demand for Deforestation per farm size & municipality
    % Recall: deforest_farm_j = [Code Area Demand]
    % IVQR
    Dem_cp = zeros(TT,size(tax_cp,2),4);
    Dem_cp(1:T(1),:,1) = deforest_farm_1cp(:,3:size(deforest_farm_1cp,2));
    Dem_cp(1:T(2),:,2) = deforest_farm_2cp(:,3:size(deforest_farm_2cp,2));
    Dem_cp(1:T(3),:,3) = deforest_farm_3cp(:,3:size(deforest_farm_3cp,2));
    Dem_cp(1:T(4),:,4) = deforest_farm_4cp(:,3:size(deforest_farm_4cp,2));

    % 2SLS
    Dem_cp_iv = zeros(TT,size(tax_cp,2),4);
    Dem_cp_iv(1:T(1),:,1) = deforest_farm_1cp_iv(:,3:size(deforest_farm_1cp_iv,2));
    Dem_cp_iv(1:T(2),:,2) = deforest_farm_2cp_iv(:,3:size(deforest_farm_2cp_iv,2));
    Dem_cp_iv(1:T(3),:,3) = deforest_farm_3cp_iv(:,3:size(deforest_farm_3cp_iv,2));
    Dem_cp_iv(1:T(4),:,4) = deforest_farm_4cp_iv(:,3:size(deforest_farm_4cp_iv,2));
    
    % Area
    % Recall: Data_farm = [Code Y Area Qa1 Qa2 Uivqr beta_ivqr(1,Uivqr)']; 
    Areacp = zeros(TT,4);
    Areacp(1:T(1),1) = Data_farm_1(:,3);
    Areacp(1:T(2),2) = Data_farm_2(:,3);
    Areacp(1:T(3),3) = Data_farm_3(:,3);
    Areacp(1:T(4),4) = Data_farm_4(:,3);

    % Coefficients from IVQR estimation
    beta_ivqr = zeros(TT,4);
    beta_ivqr(1:T(1),1) = Data_farm_1(:,7);
    beta_ivqr(1:T(2),2) = Data_farm_2(:,7);
    beta_ivqr(1:T(3),3) = Data_farm_3(:,7);
    beta_ivqr(1:T(4),4) = Data_farm_4(:,7);
   
elseif penalty == 1
    
    % Sample size for each farm size
    T(1) = size(deforest_farm_1cp_fines,1);
    T(2) = size(deforest_farm_2cp_fines,1);
    T(3) = size(deforest_farm_3cp_fines,1);
    T(4) = size(deforest_farm_4cp_fines,1);
    TT   = max(T);

    % Demand for Deforestation per farm size & municipality
    % Recall: deforest_farm_j = [Code Area Demand]
    % IVQR
    Dem_cp = zeros(TT,size(tax_cp,2),4);
    Dem_cp(1:T(1),:,1) = deforest_farm_1cp_fines(:,3:size(deforest_farm_1cp_fines,2));
    Dem_cp(1:T(2),:,2) = deforest_farm_2cp_fines(:,3:size(deforest_farm_2cp_fines,2));
    Dem_cp(1:T(3),:,3) = deforest_farm_3cp_fines(:,3:size(deforest_farm_3cp_fines,2));
    Dem_cp(1:T(4),:,4) = deforest_farm_4cp_fines(:,3:size(deforest_farm_4cp_fines,2));

    % 2SLS
    Dem_cp_iv = zeros(TT,size(tax_cp,2),4);
    Dem_cp_iv(1:T(1),:,1) = deforest_farm_1cp_iv_fines(:,3:size(deforest_farm_1cp_iv_fines,2));
    Dem_cp_iv(1:T(2),:,2) = deforest_farm_2cp_iv_fines(:,3:size(deforest_farm_2cp_iv_fines,2));
    Dem_cp_iv(1:T(3),:,3) = deforest_farm_3cp_iv_fines(:,3:size(deforest_farm_3cp_iv_fines,2));
    Dem_cp_iv(1:T(4),:,4) = deforest_farm_4cp_iv_fines(:,3:size(deforest_farm_4cp_iv_fines,2));
    
    % Area
    % Recall: Data_farm = [Code Y Area Qa1 Qa2 Uivqr beta_ivqr(1,Uivqr)']; 
    Areacp = zeros(TT,4);
    Areacp(1:T(1),1) = Data_farm_1(:,3);
    Areacp(1:T(2),2) = Data_farm_2(:,3);
    Areacp(1:T(3),3) = Data_farm_3(:,3);
    Areacp(1:T(4),4) = Data_farm_4(:,3);

    % Coefficients from IVQR estimation
    beta_ivqr = zeros(TT,4);
    beta_ivqr(1:T(1),1) = Data_farm_1(:,7);
    beta_ivqr(1:T(2),2) = Data_farm_2(:,7);
    beta_ivqr(1:T(3),3) = Data_farm_3(:,7);
    beta_ivqr(1:T(4),4) = Data_farm_4(:,7);

    
end

% Bandwidth for local approximations
h = 0.01;

% Matrices to be Used
% IVQR
tax_hat_cp = zeros(TT,4); % Local approx for tax
dem_hat_cp = zeros(TT,4); % Local approx for corresponding demand function
LS_cp      = zeros(TT,4); % Lost Surplus for each farm size and municipality

% 2SLS
tax_hat_cp_iv = zeros(TT,4); % Local approx for tax
dem_hat_cp_iv = zeros(TT,4); % Local approx for corresponding demand function
LS_cp_iv      = zeros(TT,4); % Lost Surplus for each farm size and municipality


% II. Compute Lost Surplus for each municipality and each farm size

% Fix Farm Size
for f = 1:4
                
    % Fix Municipality
    for j = 1:T(f)
        
        % 1. IVQR
        if beta_ivqr(j,f)< 0
            
            % (a) Share of Agriculture
            share = Dem_cp(j,:,f)./Areacp(j,f); 
            
            % (b) Local Approx using kernel at "prop" - the proportion
            %     aimed by the policy
            aprox = ksrmv(share',tax_cp',h,prop);
            
            % Change bandwidth if local approximation does not work
            while (isnan(aprox.f)==1)
                h = h + 0.01;
                aprox = ksrmv(share',tax_cp',h,prop);
            end
            h = 0.01;
            
            % (c) Predicted Tax and Share from policy
            tax_hat_cp(j,f) = aprox.f;
            dem_hat_cp(j,f) = prop*Areacp(j,f);
            
            % (d) Find position of the predicted tax
            m = 1;
            while (tax_cp(m) < tax_hat_cp(j,f) && m < size(tax_cp,1))
                m = m + 1;
            end
            
            % (e) Vector of the Demand Function up to the "prop" point
            if m == 1
                tax_star_cp = [0             tax_hat_cp(j,f)];
                dem_star_cp = [Dem_cp(j,1,f) dem_hat_cp(j,f)];
            else
                tax_star_cp = [tax_cp(1,1:m-1)   tax_hat_cp(j,f)];
                dem_star_cp = [Dem_cp(j,1:m-1,f) dem_hat_cp(j,f)];
            end
                                  
            % (f) Compute Lost Surplus for farm size f in muncipality j
            LS_cp(j,f) = trapz(tax_star_cp, dem_star_cp);
            
        end
        
        % 2. 2SLS
        
            % (a) Share of Agriculture
            share_iv = Dem_cp_iv(j,:,f)./Areacp(j,f); 
            
            % (b) Local Approx using kernel at "prop" - the proportion
            %     aimed by the policy
            aprox_iv = ksrmv(share_iv',tax_cp',h,prop);
            
            % Change bandwidth if local approximation does not work
            while (isnan(aprox_iv.f)==1)
                h = h + 0.01;
                aprox_iv = ksrmv(share_iv',tax_cp',h,prop);
            end
            h = 0.01;
            
            % (c) Predicted Tax and Share from policy
            tax_hat_cp_iv(j,f) = aprox_iv.f;
            dem_hat_cp_iv(j,f) = prop*Areacp(j,f);
            
            % (d) Find position of the predicted tax
            m = 1;
            while (tax_cp(m) < tax_hat_cp_iv(j,f) && m < size(tax_cp,1))
                m = m + 1;
            end
            
            % (e) Vector of the Demand Function up to the "prop" point
            if m == 1
                tax_star_cp_iv = [0                tax_hat_cp_iv(j,f)];
                dem_star_cp_iv = [Dem_cp_iv(j,1,f) dem_hat_cp_iv(j,f)];
            else
                tax_star_cp_iv = [tax_cp(1,1:m-1)      tax_hat_cp_iv(j,f)];
                dem_star_cp_iv = [Dem_cp_iv(j,1:m-1,f) dem_hat_cp_iv(j,f)];
            end
                                              
            % (f) Compute Lost Surplus for farm size f in muncipality j
            LS_cp_iv(j,f) = trapz(tax_star_cp_iv, dem_star_cp_iv);
                    
    end
    
end

% III. Compute aggregated lost surplus

% IVQR
% Total Lost Surplus
LS_80rule_cp = sum(sum(LS_cp));

% 2SLS
% Total Lost Surplus
LS_80rule_cp_iv = sum(sum(LS_cp_iv));


%% (D)  SHOW RESULTS

% I. Show Lost Surplus Results

if results_est == 1
    
    disp('----------------------------------------------------------------------------')
    disp(' ')
    disp('LOST SURPLUSES (U$ million) - IVQR')
    table(LS_80rule_cp/10^6, LS_uniform_tax_cp/10^6, tax_revenue_cp/10^6, area_below_demand_cp/10^6,...
        'VariableNames',{'Rule_80Percent_LostSurplus','UniformTax_LostSurplus',...
        'UniformTax_TaxRevenues','UniformTax_AreaBelowDemand'},...
        'RowNames',{'Estimates'})
    
    % II. Show Uniform Tax
    if penalty == 0     % Actual Monitoring       
           
        disp('----------------------------------------------------------------------------')
        disp(' ')
        disp('VALUE OF THE UNIFORM TAX (U$/ha) & CORRESPONDING LOST SURPLUSES (U$ million) - IVQR & 2SLS')
        table([tax_ivqr tax_iv]',...
            [LS_uniform_tax_cp/10^6 LS_uniform_tax_cp_iv/10^6]',...    
            'VariableNames',{'Value_UniformTax','UniformTax_Lost_Surplus'},...
            'RowNames',{'IVQR','2SLS'})
        disp('----------------------------------------------------------------------------')
        
    elseif penalty == 1 % Set fines at pre-2004 levels
        
        disp('----------------------------------------------------------------------------')
        disp(' ')
        disp('VALUE OF THE UNIFORM TAX (U$/ha) & CORRESPONDING LOST SURPLUS (U$ million) - IVQR & 2SLS')
        table(tax_ivqr,...
            LS_uniform_tax_cp/10^6,...    
            'VariableNames',{'Value_UniformTax','UniformTax_Lost_Surplus'},...
            'RowNames',{'IVQR'})
        disp('----------------------------------------------------------------------------')
        
    end
    
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
